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Numerical analysis of co-firing pulverized coal and biomass in a vertical cylindrical laboratory furnace is 
explored. The ratio of coal and biomass in the fuel was 80:20 by mass for all cases. The mathematical 
model of combustion in the furnace was established by describing physical phenomena such as turbulent 
flow, heat and mass transfer, devolatilization and combustion. 

A 3D-model of combustion in a laboratory furnace was created using the CFD software FLUENT. The 
shape of the biomass particles was estimated as cylindrical and was accounted for in the calculation of 
particle trajectories via a custom-developed model. Experimental measurements were conducted on a 
20 kW laboratory furnace with controllable wall temperature. The temperature varied in the range from 
1233 K to 1823 K, depending on the case. Excess air for combustion was set at 10% or 20%, depending on 
the case. The developed model shows better agreement with the experimental data than the existing 
models, which estimate particles as spheres. Analysis of the results shows that the influence of the 
particle size increases with the size of the particle. Also, the geometry of the cylindrical particles strongly 
influences the beginning and the intensity of the devolatilization process and subsequently the com¬ 
bustion process. 

© 2014 Elsevier Ltd. All rights reserved. 


1. Introduction 

Mechanisms of coal combustion in industrial furnaces are well 
understood today. Many mathematical and numerical models are 
developed on the basis of these insights that successfully describe 
the observed physical processes such as in Refs. [1—3]. However, it 
is unrealistic to expect that these models are applicable in the 
modelling of biomass combustion since biomass somewhat differs 
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from coal in its physical and chemical properties. The variety of 
biomass particle shapes and sizes plays an important role in its 
combustion process as different particle surface to volume ratios 
directly affect the rate of the formation of volatiles and oxidation 
during combustion. The pulverized biomass particle can usually be 
described as a cylinder or disc. Both of these geometric shapes have 
a much higher surface to volume ratio than the spheres which are 
commonly used to describe coal dust particles. As a consequence, it 
is necessary to pay special attention to the influence of biomass 
particle shape and size when creating a biomass combustion 
model. 

Most existing studies are applicable only to the individual par¬ 
ticles or particles of a spherical shape and do not involve com¬ 
bustion. Only recently have a few studies appeared which include 














I. Bonefacic et al. / Applied Thermal Engineering 78 (2015) 74-81 


75 


the combustion processes. In Ref. [4] the effect of temperature and 
particle size on biomass pyrolysis and combustion of char residue 
was studied. One step further was made in Ref. [5] and combustion 
of individual biomass particles of different shapes and of the same 
volume was modelled, which was described as disks, cylinders, or 
as nearly spherical. The combustion of pulverized biomass or the 
co-firing of coal and biomass in furnaces is generally modelled with 
spherical particles where, in some cases, unsphericity is taken into 
account only through the form factors. Thus, in Backreedy et al. [6] 
or in more recent papers [7,8], the drag coefficient for the particle 
flow model by Haider and Levenspiel [9] was used to model pul¬ 
verized coal and biomass combustion. 

The aim of this paper is to set a suitable mathematical model 
that would more accurately describe the physical process of co¬ 
firing pulverized coal and biomass. Special attention is given to 
the modelling of irregularly shaped biomass particles which are 
approximated by cylinders. The results obtained by numerical 
simulation are compared with the measurements performed on a 
20 kW entrained flow reactor designed for co-firing pulverized coal 
and biomass, installed at the Faculty of Mechanical Engineering, 
University of Sarajevo. 

2. Mathematical model 

The commercial CFD code FLUENT 6.3 [10] was used with a 
number of parameters appropriate for coal combustion. CFD 
models for coal combustion contain a number of sub-models for 
turbulent flow, heat transfer and combustion of the coal. Many of 
these processes are common to the combustion of biomass but 
some aspects can be significantly different, particularly the particle 
shape and size, content and the nature of the volatiles and the ki¬ 
netics of the devolatilisation process. 

In this work, release of volatiles from the coal particles is 
described by the Two Competing Rates model [11 ], while the release 
of volatiles from biomass particles is described by the somewhat 
simpler Single Kinetic Rate model [12] due to the lack of necessary 
data for biomass such as the volatile release rate and the exact 
composition of the released volatiles. 

Combustion of volatiles is modelled using the EDM/Finite-Rate 
model [10]. The basic chemical reactions for combustion of coal and 
biomass volatiles have been formulated with the assumption that 
the resulting complex hydrocarbons can be described by the gen¬ 
eral formula CxHyO z S m N n . It was also necessary to include the ra¬ 
diation sub-model, as a significant portion of the heat released 
during combustion is transferred by radiation. Due to the specific 
geometry of the combustion chamber, this is done by the Discrete 
Ordinates model [13,14] which in combination with the weighted 
sum of gray gases model (WSGGM), offers the possibility of 
modelling the effects of non-gray gases as has already been proven 
in Ref. [3]. In addition to the processes taking place in the gas 
(continuous) phase, it was necessary to take into account the pro¬ 
cesses on the boundary between solid (discrete) and the gas phase. 
FLUENT contains a comprehensive model based on the Lagrangian 
approach to discrete phase combustion. User defined subroutines 
were implemented to approximate biomass cylindrical particles. 
The developed model, in addition to existing drag forces on coal 
spherical particles in turbulent flow, calculates the drag forces and 
lift forces on cylindrical particles of biomass. Particle combustion 
was modelled using the Kinetics/Diffusion Rate-Limited model 
[15,16] that takes into account the speed of oxidation on the surface 
of the particles and the rate of the diffusion of oxygen to the surface 
of the particles and takes the slower of the two as relevant. In post¬ 
processing, based on the calculated velocity and temperature field 
and concentration of the species, the pollutant concentrations NO 
and SO 2 are calculated. 


2.1. Cylindrical particle model 


Particle models available in FLUENT calculate drag force on the 
particle independently of the particle orientation; though lift force is 
not included, which is somewhat acceptable for modelling spherical 
particles such as of pulverized coal. Biomass particles, however, are 
of irregular geometry and those factors need to be included. It was, 
therefore, necessary to create a model of biomass particle trajectory 
calculation in turbulent flow that approximates its shape as a cyl¬ 
inder or a disc, including the lift forces acting upon it. Coal particles 
were, at the same time, approximated as spheres with models 
available in FLUENT where drag force on the particle is defined as: 


Fd = ?C DP A ef i\~u - Ttp\(u - Up). 


(1) 


where Cd is the drag coefficient and A e p\ is the effective area rep¬ 
resenting the particle projected area in the direction of the drag 
force. For the spherical particle, it is equal to the circular area of the 
diameter d p . 

According to [17], where a large number of available correlations 
for drag coefficient calculation is processed, the model by Ganser 
[18], which takes into account the actual shape of the particle, has 
proven to be the most accurate. The effective area for a cylinder of 
length L and a diameter d p is: 


d 2 /- 

A e f ^ = 7ry cos 2 (4L/d p 7r) 2 sin 2 a. (2) 

It is dependant on the angle between the relative particle ve¬ 
locity (Tt - Ti p ) and its axis ~z' (Fig. 1). 

The drag coefficient for the cylindrical particle defined by 
Ganser [18] is: 
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Fig. 2. Schematic layout of the entrained flow reactor, [20]. 


where the Reynolds number, Re s /, is calculated for a sphere of equal 
volume. Coefficients JCi and I< 2 are defined as: 

Kl = (wf + W*) ’ l0 SK 2 = l-8148(—log 0) 05743 , (4) 

where 0 is the shape factor defined as the ratio of the surface area of 
a sphere having the same volume as the particle to the actual 
surface area of the particle, d s /is the diameter of a sphere having the 
same volume as the particle and d n is the diameter of a circle having 
an area equal to the effective area A e p. Beside drag force, it is also 
necessary to include lift force on the particle due to its cylindrical 
shape (Fig. 1). 


| Primary air+fuel 



Lift force F L is perpendicular to the particle relative velocity 
(u - Tip) and with appropriate simplifications to the drag force 
F D . Also, lift force lies in a plane defined by the particle axis ~z and 
relative velocity (it - Tt p ), therefore it is defined as: 

Fl = ^C LP A ef 2 —A —=r-p2- [/ x (u - Up)] x (u - u p ), 

(5) 

where A e f i2 is the particle projected area in the direction of the lift 
force: 

d 2 /- 

A e f 2 = 7ry sin 2 a + (4L/d p n) 2 cos 2 a. (6) 

The intensity of the lift force can be determined from the 
expression [19],: 

\fj\ I 9 I 

4^- = sin a*cos a (7) 

\Fd\ 1 1 



Fig. 3. Furnace and burner geometry discretization. 


Fig. 4. Macrophotography of biomass particles in tested fuel. 
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Table 1 

Fuel particle size distribution, rest on the sieve in %wt. 


Fuel 

>500 pm 

200-500 pm 

100-200 pm 

45-100 pm 

<45 pm 

Coal [%] 

0 

18.2 

28.0 

17.4 

36.4 

Biomass [%] 

0 

67.5 

26.7 

4.2 

1.6 


When the angle of attack a is not zero, the centre of pressure on 
the particle surface in motion does not coincide with the centre of 
its mass; therefore, momentum occurs that tries to rotate the 
particle. For simplicity, this momentum is ignored in this paper but 
rotation of the particle is taken into account with stochastic 
tracking in trajectory calculation and a large number of trajectories. 

2.2. Momentum exchange 

Momentum transfer from the continuous phase to a particle 
when passing through the control volume is calculated using the 
expression: 

~F = T D + ~Fg + ~F PC + ~F L , ( 8 ) 

where forces acting on the particle are: drag force, force of gravity, 
force due to pressure gradient and lift force. In the case of spherical 
particles, lift force is zero. 

3. Model validation 

The presented mathematical model of pulverized solid fuel 
combustion involving a developed sub-model for cylindrical 
biomass particle trajectory calculation was verified by comparison 
with experimental data. Experimental testing was conducted on a 
AUCR reactor ( Automatically controlled tubular reactor ) for com¬ 
bustion of pulverized solid fuels with a nominal power of 20 kW, 
designed and installed at the Faculty of Mechanical Engineering, 
University of Sarajevo, Fig. 2. It is a flexible combustion chamber in 
terms of setting the desired combustion temperature, as well as 
changes in the distribution of combustion air. The reactor is used 
for the examination of combustion characteristics of different types 
of coal and coal-biomass blends. It consists of a central vertical tube 
made of porous alumino-silicate ceramic in which combustion 
takes place. The length of the ceramic tube is 2900 mm, the 
diameter is 230/200 mm and it consists of 15 interconnected parts. 
Furnace wall temperatures up to 1833 K are possible and are 
maintained by SiC heaters with a maximum power of 70 kW, 
divided into four zones. The wall temperature of the furnace is 
measured by WIKA ptl0% Rh -1600 °C thermocouples and is 
regulated by Siemens Simatic PLC devices in each of the four heated 
zones. After the blower, the combustion air is divided on the pri¬ 
mary air, which also acts as fuel carrier, secondary air and tertiary 
air and leads to a burner located at the top of the furnace. 

Based on the furnace geometry, only one quarter of the domain 
was modelled and discretized with 120,000 control volumes. The 


Table 3 

Fuel proximate and ultimate analysis. 

Brown coal from “Kakanj” Biomass (spruce sawdust) 

Prox. analysis [%]-ar 
Moisture 
Ash 

Volatiles 
Char 


Ult. analysis [%]-ar 


C 

35.11 

38.91 

H 

2.78 

7.36 

O 

7.03 

41.61 

S 

2.28 

0.33 

N 

0.96 

0.33 

Heating val. 

. [kj/kg]-ar 


Gross 

13,490 

17,386 

Net 

12,657 

15,612 


entire furnace burner was modelled separately with 470,000 vol¬ 
umes to achieve a more accurate inlet flow profile, Fig. 3. 

When calculating the trajectories of the discrete phase, the 
shape of the coal particle is simplified by the spherical shape while 
biomass particles are assumed to be of cylindrical shape and are 
modelled with a custom-developed model. The shape of a cylinder 
with a length to diameter ratio, L/D = 10, was used as a simplifi¬ 
cation based on the actual shape of the biomass particles visible 
from the fuel sample macrophotography, Fig. 4. 

There are four input streams of particles due to their size ac¬ 
cording to the values in Table 1 : Stream of coal particles 40 pm in 
diameter; Stream of coal particles 150 pm in diameter; Stream of 
biomass particles 150 pm in length; Stream of biomass particles 
300 pm in length. Each stream is defined with 50 trajectories of 
representative particles released from 20 points at the inlet 
boundary making a total of 1000 trajectories for each of the defined 
particle sizes. 

The devolatilisation stage begins instantly after entering the 
domain and is characterized by the release of light gases (CO, CO 2 , 
FI 2 O, H 2 , etc.) as well as light and heavy hydrocarbons from coal and 
biomass as shown in the following reactions: 

Biomass -► Biomass volatiles + Char 

Coal -► Coal volatiles + Char 

The values of pre-exponential factors and the activation en¬ 
ergies of coal (Two Competing Rates Model ) and biomass ( Single 
Kinetic Rate Model ) devolatilisation are given in Table 2. The 
apparent kinetic parameters for biomass particles of 150 pm and 
300 pm were extrapolated from the available values of character¬ 
istic particle sizes [21]. 

Volatiles released from the fuel are replaced with the molecule 
of a pseudovolatile, i.e. a molecule of a complex hydrocarbon of the 
general shape C x H y O z S m N n . It is determined from proximate and 
ultimate analysis of the fuel in Table 3. 


11.30 

11.20 

41.43 

0.26 

25.88 

75.48 

21.39 

13.06 


Table 2 

Devolatilisation model parameters. 


Biomass 

Source 

Reaction rates 

Coal 

Source 

Reaction rates 

(300 pm) 



Ai = 2.0 -10 5 s —1 



Ai = 742 s —1 

£1 = 2.83-10 7 J/kmol 

(150 pm): 

[16] 

K = Aie-^ 1 / i?r p) 

£1 = 1.046-10 8 J/kmol 
a t = 0.3 

A 2 = 1.3-10 4 s —1 

[5] 

3?! = Aie~^ /RT p^ 

Ai = 19,427 s" 1 
£1 = 3.94-10 7 J/kmol 

[16] 

K = Aie-^ 1 / i?r p) 

£2 = 1.674-10 7 J/kmol 
a 2 = 1.0 

[5] 

=A 2 e-( E2 / rt p) 
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Fig. 5. O 2 and CO 2 content in dry flue gases at A = 1.2. 


Finally, the combustion of the mixture of pseudovolatiles orig¬ 
inating from biomass and coal is defined with the following three 
reactions: 

C 2 . 851 H 9 . 674 O 3 . 446 S 0 . 014 N 0.032 +2.457 02—>2.207 CO 
+ 0.644 C0 2 + 4.837 H 2 0 + 0.014 S0 2 + 0.016 N 2 (9) 

Q. 267 H 10 . 304 O 1 . 641 S 0 . 266 N 0.256 + 4,636 O 2 ^3.303 CO 
+ 0.964 C0 2 + 5.152 H 2 0 + 0.266 S0 2 + 0.128 N 2 (10) 

CO + 0.5 0 2 —>C0 2 (11) 

Combustion in the gas phase is modelled with Finite Rate/Eddy 
Dissipation model, while combustion of char residue is represented 
by solid carbon, C( S ), which occurs on the particle surface, 

C (s) +0.5O 2 ^CO, (12) 

and is modelled with the Kinetics/Diffusion-Limited Rate model. 

4. Results and discussions 

Comparison of the results obtained by the combustion model 
with the measured values on the laboratory combustion chamber 
for all test regimes are given in Figs. 5-8. 

In the legends in Figs. 6-8, the “Cylinder model” represents the 
results obtained with the combustion model, which include the 
sub-model of the cylinder-shaped biomass particles, the “Fluent 




Fig. 6. CO content in dry flue gases at different excess air ratios. 




Fig. 7. NO x emissions in dry flue gases at different excess air ratios. 
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Fig. 8. S0 2 emissions in dry flue gases at different excess air ratios. 


model 1 ” represents the model of non-spherical particles, with the 
shape factor 0 = 0.6 available in Fluent [9], and the “Fluent model 
2 ” represents the model of spherical particles also available in 
Fluent [22]. 

From the diagram in Fig. 6 , improvement in the predicted 
values of CO in the flue gases at lower test temperatures 
(T = 1413 K—1673 I<) is visible using the developed model of cy¬ 
lindrical particles as compared to the standard models. Only at the 
highest test temperature (T = 1823 K), where heating rates of the 
particles are the highest and thus the process of volatile release 
and the combustion reactions are more rapid, the standard model 
of spherical particles is closer to predicted values of the developed 
model. Significant deviations from predicted values of O 2 , CO 2 and 
CO concentrations in the flue gases can be observed at the lowest 
test temperature of 1233 K (Figs. 5 and 6 ). Such a result can be 
partly attributed to insufficiently accurate models that describe 
processes of devolatilization (two competing rates and single ki¬ 
netic rate models) and combustion (EDM/Finite rate model). The 
aforementioned models imply fast chemical reactions whereas at 
lower temperatures these processes take place more slowly and 
are more dependent on the type and composition of the fuel. To 
make the model of combustion more precise, at lower tempera¬ 
tures it should apply some of the more complex models of devo¬ 
latilization and combustion. Such models require detailed 
knowledge of the structure of the fuel and the exact composition 
of volatiles, primarily biomass, but such experimental data is very 
hard to come by. 

NO and SO 2 emissions were measured and modelled in ppm 
converted to mg/m„ and normalized at 6%02 in dry combustion 
gases. When modelling the emission of nitric oxide thermal NO and 
fuel, NO mechanisms were taken into account. Fig. 7 shows the 
calculated values of the concentration of nitrogen monoxide for all 
test temperatures. In test cases with an excess air ratio of 1.2 and 
compared with the measured values, increase in the concentration 
of NO at the highest process temperature (1823 I<) can be observed. 
It is a result of thermal NO which is substantially formed by the 
Zeldovich mechanism above 1673 K. Also, a drop in the concen¬ 
tration of nitrogen monoxide can be noticed by comparing the di¬ 
agrams in Fig. 7 as a result of excess air reduction. 

Predicted NO concentrations are slightly lower than measured, 
however, and the trend of changes with temperature is consistent 
with all compared models which can be attributed to the known 
temperature field. It is possible to perceive a somewhat better 
agreement with measured concentrations of NO achieved using a 
combustion model that takes into account the cylindrical shape of 
biomass particles. The values given by the spherical model are 
significantly lower than the measured values, probably as a result 


of the combustion of fuel particles which are predicted as insuf¬ 
ficiently scattered. In such conditions, a fuel-rich zone is created in 
which the fuel nitrogen has a tendency to form molecular N 2 . 

Overall, results obtained by the presented combustion model 
coincide with the measured values in the laboratory furnace for 
temperatures above 1400 K, except for SO 2 emissions in flue gases. 
Measured emissions of sulphur dioxide (Fig. 8 ) were significantly 
less than the values predicted by the model. A possible reason for 
the excessive proportion of predicted SO 2 is the desulfation process 
which is not taken into account by the developed model. The cause 
of the desulfation process in the combustion chamber is the very 
high ash content in brown coal “Kakanj” and, at the same time, the 
high content of CaO in the ash of about 12%. Fig. 9 shows the 
relationship between the sulphur capture rate Sb, defined as the 
portion of SO 2 that reacts with basic oxides from the ashes (mainly 
CaO) and the temperature in the combustion chamber for tested 
fuels. 

The nearly identical predicted values of SO 2 emissions obtained 
by all three discrete phase models are expected given that all 
sulphur in the tested fuel originates from coal particles (as can be 
seen from fuel ultimate analysis in Table 3), which are approxi¬ 
mated by spheres in all three models. 

The calculated trajectories of the biomass particles defined as 
cylinders with a length of 150 pm and 300 pm and of a sphere of 
equal volume can be seen in Fig. 10. It can be seen that biomass 
particles are significantly more scattered when modelled as cylin¬ 
ders, which is especially pronounced with larger particles. Differ¬ 
ences in calculated particle trajectories lead to differences in the 
rate of volatile release and its concentration in the combustion 
chamber (Fig. 11 ) and as a result, the calculated composition of the 
flue gases. 



t, °C 


Fig. 9. Sulphur capture rate Sb as a function of process temperature for the “Kakanj” 
coal (K)- spruce sawdust (S), co-firing at different ratios (K:S) at A = 1.2, [23]. 
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Cylinder (Z,=150 pm) 



Sphere 



Cylinder (1=300 pm) 



Sphere 



Fig. 10. Comparison of predicted trajectories for biomass particles modelled as cylinders or spheres of equal volume for the case T = 1673 K and A = 1.2. 


5. Conclusion 

The developed model describing biomass particles as cylinders 
proved to be more accurate for the prediction of carbon monoxide 
and nitrogen monoxide concentrations in flue gases than existing 
models that describe particles as spheres. Given that the test is 
conducted using a mixture of coal and biomass in a weight ratio of 
80:20 it is reasonable to assume that this difference would be even 
greater if a fuel mixture with a higher proportion of biomass was 
used for testing. 

The model was subsequently used to analyze the influence of 
the geometry of the particles in the combustion process. It revealed 
that geometry has a significant impact on the starting point and the 
rate of devolatilisation as well as the burning of coke residue. The 
higher length to diameter ratio of the cylinder that describes the 


particle leads to increased scattering of fuel particles after entering 
the domain. This extends the devolatilization zone; the process of 
devolatilization is accelerated and the mixing of volatiles with the 
air is better. Fuel particles modelled by sphere are much less scat¬ 
tered and the concentration of particles is increased along the 
chamber axis. The result is combustion of volatiles in the fuel-rich 
zone, which affects local temperatures, concentrations of com¬ 
bustion products and afterwards the combustion of the coke 
residue. 
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Fig. 11. Predicted devolatilization rates from biomass particles modelled as cylinders 
or spheres of equal volume for the case T = 1673 K and A = 1.2. 
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